#################################################################
##                        Preliminaries                        ##
#################################################################

# Load useful packages  -----------------------------------------------------------

# install.packages("pacman") # Uncomment if `pacman` package is not already
# installed
library(pacman)

p_load(
tidyverse,
janitor,
broom,
countrycode,

# Modeling
lmtest,
plm,
tscsdep,
spdep,
spatialreg,
urca,
dynamac,
fixest,
bayestestR,

# Graphs and Figures
patchwork,
marginaleffects,

# Tables
texreg,
modelsummary,
kableExtra,

# Project managements
here
)

# Tell R where you are
i_am("01_script_master_women_exchange_replication.R")

# Load useful functions  ---------------------------------------------


# Create a function that displays coefficient estimates with PCSEs
show_pcse <- function (x) {
  #'  Show Beck and Katz (1995) panel-corrected standard errors standard errors
  #'  for time-series cross-section linear models. Note that this code will 
  #'  *only* work with models estimated with the `plm` package.
  
  #' @param x:  A `plm` linear model object.
  #'    
  #' @return The `summary.plm` method for plm objects but with the Beck and 
  #' Katz (1005) standard errors and associated t-values and p-values.
  #'
  #' @example show_pcse(x)
  
  summary(x, vcovBK(x, type = "HC1", cluster = "time"))
}

# I use `texreg` to make the tables. When using "robust" standard errors, 
# `texreg` requires functions to extract the correct standard errors and
# p-values. The two functions extract the relevant items.
extract_pcse <- function(x) {
  #' Extract the Beck and Katz (1995) standard errors for use in `texreg`
  #' tables. Note that this code will *only* work with model estimated with the
  #' `plm` package.
  
  #' @argument x: A plm linear model object.
  
  #' @returns When used inside the `texreg` function with the standard error
  #'   specified with the `override.se` list, this function allows the user to
  #'   substitute OLS standard errors with those recommended by Beck and Katz
  #'   (1995).
  #'
  #' @example texreg(x, override.se = extract_pcse(x))

  coeftest(x, vcov. = vcovBK(x, type = "HC1", cluster = "time"))[, 2]
  
}

extract_pcse_sig <- function(x) {
  #' Extract the p-value associated with the Beck and Katz (1995) standard
  #' errors for use in `texreg` tables Note that this code will *only* work
  #' with model estimated with the `plm` package.
  
  #' @argument x: A plm linear model object.
  
  #' @returns When used inside the `texreg` function with the p-values specified
  #'   with the `override.pvalue` list, this function allows the user to
  #'   substitute the "stars" for OLS standard errors with those commensurate
  #'   with the Beck and Katz (1995) standard errors.
  #'   
  #'   @example texreg(x, override.se = extract_pcse(x), override.pvalues = extract_pcse_sig(x))
  
  coeftest(x, vcov. = vcovBK(x, type = "HC1", cluster = "time"))[, 4]
}

# One model in the paper uses the `dynardl` package, but `texreg` doesn't
# support those objects out of the box. This can be fixed by adding a custom
# `extract` function that grabs all the relevant information.
extract.dynardl <- function(model, 
                            include.rsquared = FALSE, 
                            include.adj.rsquared = FALSE, 
                            include.nobs = TRUE) {
  s <- summary(model)
  names <- rownames(s$coefficients)
  co <- s$coefficients[, 1]
  se <- s$coefficients[, 2]
  pval <- s$coefficients[, 4]
  
  gof <- numeric()
  gof.names <- character()
  gof.decimal <- logical()
  
  if (include.rsquared == TRUE) {
    rs <- s$r.squared
    gof <- c(gof, rs)
    gof.names <- c(gof.names, "R$^2$")
    gof.decimal <- c(gof.decimal, TRUE)
  }
  
  if (include.adj.rsquared == TRUE) {
    adj <- s$adj.rsquared
    gof <- c(gof, adj)
    gof.names <- c(gof.names, "Adj.\\ R$^2$")
    gof.decimal <- c(gof.decimal, TRUE)
  }
  
  if (include.nobs == TRUE) {
    n <- nrow(model$model$model)
    gof <- c(gof, n)
    gof.names <- c(gof.names, "Num.\\ obs.")
    gof.decimal <- c(gof.decimal, TRUE)
  }

    tr <- createTexreg(
    coef.names = names,
    coef = co,
    se = se,
    pvalues = pval,
    gof.names = gof.names,
    gof = gof
    
  )
  return(tr)
}

# Now, we register the function as a method for the generic extract() function
setMethod("extract",
          signature = className("dynardl", "dynardl"),
          definition = extract.dynardl)

# Load useful settings for tables -----------------------------------------

# Make list of variable names for `texreg` tables
list_var_names_texreg <- list(
  'flfplag'                           = 'LDV$_{(t-1)}$',
  'flfp_mlfplag'                      = 'LDV$_{(t-1)}$',
  'dflfplag'                          = '$\\Delta$LDV$_{(t-1)}$',
  'dflfp_mlfplag'                     = '$\\Delta$LDV$_{(t-1)}$',
  'log_wdi_overvaluedlag'             = 'logOvervalued$_{(t-1)}$',
  'dlog_wdi_overvalued'               = '$\\Delta$logOvervalued$_t$', 
  'log_gdpcaplag'                     = 'GDP$_{(t-1)}$',  
  'dlog_gdpcap'                       = '$\\Delta$GDP$_t$',  
  'sqlog_gdpcaplag'                   = 'GDP$^{2}_{(t-1)}$',  
  'dsqlog_gdpcap'                     = '$\\Delta$GDP$^{2}_t$',  
  'log_oilgascaplag'                  = 'Resource Rents$_{(t-1)}$',  
  'dlog_oilgascap'                    = '$\\Delta$Resource Rents$_t$',
  'wdi_fertilitylag'                  = 'Fertility$_{(t-1)}$',
  'dwdi_fertility'                    = '$\\Delta$Fertility$_t$',
  'vdem_elec_demlag'                  = 'Regime Type$_{(t-1)}$', 
  'dvdem_elec_dem'                    = '$\\Delta$Regime Type$_t$',
  'as_factor_irr_rate_regimelag2'     = 'Exchange Rate: Narrow Crawling$_t$',
  'as_factor_irr_rate_regimelag3'     = 'Exchange Rate: Wide Crawling$_t$',
  'as_factor_irr_rate_regimelag4'     = 'Exchange Rate: Freely Floating$_t$',
  'as_factor_irr_rate_regimelag5'     = 'Exchange Rate: Freely Falling$_t$',
  'as_factor_irr_rate_regimelag6'     = 'Exchange Rate: Dual Market$_t$',
  'rho'                               = 'Spatial Lag: $\\Delta$Dependent Variable$_t$', 
  '$\\rho$'                           = 'Spatial Lag: $\\Delta$Dependent Variable$_t$'
)


# Load useful settings for figures ----------------------------------------

# Figures
theme_set(theme_minimal(base_size = 10))
theme_update(text = element_text(family = "sans"))


# Load baseline data --------------------------------------------------------

# TSCS data frame spanning 1990--2015 inclusive
df_tscs <-
  read_csv("05_data_women_exchange_tscs_replication.csv") %>% 
  group_by(cowcode) %>% 
  # Create lags and differenced variables
  mutate(dflfplag = dplyr::lag(dflfp, 1),
         flfplag2 = dplyr::lag(flfp, 2),
         dflfp_mlfplag = dplyr::lag(dflfp_mlfp, 1)) %>% 
  ungroup() 


# The spatial regressions need alternative names for the some countries to 
# accurately create k-nearest neighbors. Fix those names here.
df_tscs  <- mutate(df_tscs, country_space = country)

df_tscs$country_space <- recode_factor(df_tscs$country_space, "Bahamas, The"="Bahamas")
df_tscs$country_space <- recode_factor(df_tscs$country_space, "Belarus"="Belarus (Byelorussia)")
df_tscs$country_space <- recode_factor(df_tscs$country_space, "Bosnia and Herzegovina"="Bosnia-Herzegovina")
df_tscs$country_space <- recode_factor(df_tscs$country_space, "Brunei Darussalam"="Brunei")
df_tscs$country_space <- recode_factor(df_tscs$country_space, "Burkina Faso"="Burkina Faso (Upper Volta)")
df_tscs$country_space <- recode_factor(df_tscs$country_space, "Cabo Verde"="Cape Verde")
df_tscs$country_space <- recode_factor(df_tscs$country_space, "Cambodia"="Cambodia (Kampuchea)")
df_tscs$country_space <- recode_factor(df_tscs$country_space, "Congo, Dem. Rep."="Congo, Democratic Republic of (Zaire)")
df_tscs$country_space <- recode_factor(df_tscs$country_space, "Congo, Rep."="Congo")
df_tscs$country_space <- recode_factor(df_tscs$country_space, "Cote d'Ivoire"="Cote D'Ivoire")
df_tscs$country_space <- recode_factor(df_tscs$country_space, "Egypt, Arab Rep."="Egypt")
df_tscs$country_space <- recode_factor(df_tscs$country_space, "Ethiopia -pre 1993"="Ethiopia")
df_tscs$country_space <- recode_factor(df_tscs$country_space, "Ethiopia 1993-"="Ethiopia")
df_tscs$country_space <- recode_factor(df_tscs$country_space, "Gambia, The"="Gambia")
df_tscs$country_space <- recode_factor(df_tscs$country_space, "Germany"="German Federal Republic")
df_tscs$country_space <- recode_factor(df_tscs$country_space, "Iran, Islamic Rep."="Iran (Persia)")
df_tscs$country_space <- recode_factor(df_tscs$country_space, "Italy"="Italy/Sardinia")
df_tscs$country_space <- recode_factor(df_tscs$country_space, "Korea, Rep."="Korea, Republic of")
df_tscs$country_space <- recode_factor(df_tscs$country_space, "Lao PDR"="Laos")
df_tscs$country_space <- recode_factor(df_tscs$country_space, "Macedonia, FYR"="Macedonia (FYROM/North Macedonia)")
df_tscs$country_space <- recode_factor(df_tscs$country_space, "Madagascar"="Madagascar (Malagasy)")
df_tscs$country_space <- recode_factor(df_tscs$country_space, "Myanmar"="Myanmar (Burma)")
df_tscs$country_space <- recode_factor(df_tscs$country_space, "Pakistan-post-1972"="Pakistan")
df_tscs$country_space <- recode_factor(df_tscs$country_space, "Pakistan-pre-1972"="Pakistan")
df_tscs$country_space <- recode_factor(df_tscs$country_space, "Romania"="Rumania")
df_tscs$country_space <- recode_factor(df_tscs$country_space, "Russian Federation"="Russia (Soviet Union)")
df_tscs$country_space <- recode_factor(df_tscs$country_space, "Slovak Republic"="Slovakia")
df_tscs$country_space <- recode_factor(df_tscs$country_space, "Sri Lanka"="Sri Lanka (Ceylon)")
df_tscs$country_space <- recode_factor(df_tscs$country_space, "Eswatini"="Swaziland (Eswatini)")
df_tscs$country_space <- recode_factor(df_tscs$country_space, "Suriname"="Surinam")
df_tscs$country_space <- recode_factor(df_tscs$country_space, "Syrian Arab Republic"="Syria")
df_tscs$country_space <- recode_factor(df_tscs$country_space, "Tanzania"="Tanzania (Tanganyika)")
df_tscs$country_space <- recode_factor(df_tscs$country_space, "Turkey"="Turkey (Ottoman Empire)")
df_tscs$country_space <- recode_factor(df_tscs$country_space, "United States"="United States of America")
df_tscs$country_space <- recode_factor(df_tscs$country_space, "Venezuela, RB"="Venezuela")
df_tscs$country_space <- recode_factor(df_tscs$country_space, "Vietnam"="Vietnam, Democratic Republic of")
df_tscs$country_space <- recode_factor(df_tscs$country_space, "Yemen, Rep."="Yemen (Arab Republic of Yemen)")
df_tscs$country_space <- recode_factor(df_tscs$country_space, "Zimbabwe"="Zimbabwe (Rhodesia)")



# Load ILO industry employment data -----------------------------------------------------------

source("02_script_make_ILO_data.R")


#################################################################
##                 Table 1: Summary Statistics                 ##
#################################################################


# Table 1: Summary Statistics ---------------------------------------------

# Grab the data in the baseline models to extract for summary statistics
df_sum_stats <-
  plm(
    dflfp ~ flfplag + plm::lag(dflfp, 1) + 
      log_wdi_overvaluedlag + dlog_wdi_overvalued + 
      log_gdpcaplag + dlog_gdpcap + sqlog_gdpcaplag + 
      dsqlog_gdpcap + log_oilgascaplag + dlog_oilgascap - country - year,
    data = subset(df_tscs, oecd == 0 & country != "Zimbabwe"),
    model = "within",
    effect = "individual",
    index = c("cowcode", "year")
  )$model 

# Order the variables
df_sum_stats <- df_sum_stats |>
  select(!contains("plm")) |>
  rename(
    "D.FLFP" = dflfp,
    "L.FLFP" = flfplag,
    "L.REER Overvaluation" = log_wdi_overvaluedlag,
    "D.REER Overvaluation" = dlog_wdi_overvalued,
    "L.Per Capita GDP" = log_gdpcaplag,
    "D.Per Capita GDP" = dlog_gdpcap,
    "L.Squared Per Capita GDP" = sqlog_gdpcaplag,
    "D.Squared Per Capita GDP" = dsqlog_gdpcap,
    "L.Resource Rents" = log_oilgascaplag,
    "D.Resource Rents" = dlog_oilgascap
  )

# View the table
datasummary(
  All(df_sum_stats) ~ N + Mean + SD + Min + Median + Max,
  data = df_sum_stats,
  output = "markdown",
  title = "Summary Statistics",
  notes = "L.X and D.X refer to the lagged and differenced values of X."
) 


##################################################################
##                   Table 2: Baseline Models                   ##
##################################################################


# Table 2: FLFP ~ REER  + Country Effects ---------------------------------------------------

m_flfp_base <- plm(
  dflfp ~ flfplag + dflfplag +
    log_wdi_overvaluedlag + dlog_wdi_overvalued +
    log_gdpcaplag + dlog_gdpcap +
    sqlog_gdpcaplag + dsqlog_gdpcap +
    log_oilgascaplag + dlog_oilgascap,
  data = subset(df_tscs,
                oecd == 0 &
                  country != "Zimbabwe"),
  model = "within",
  effect = "individual",
  index = c("cowcode", "year")
)

# View results with PCSEs
show_pcse(m_flfp_base)

# Sanity check that my shorthand function returns the same estimates as the
# coeftest function
coeftest(m_flfp_base,
         vcov. = vcovBK(m_flfp_base, type = "HC1", cluster = "time"))

# Table 2: FLFP ~ REER  + Extended Controls -------------------------------

# Create a dataset to estimate the model
df <- df_tscs %>% 
  filter(
    oecd == 0 &
      country_space != "Zimbabwe (Rhodesia)" &
      country_space != "Tonga" # Causes the nearest neighbor function to balk
  )

# Estimate the baseline model using `lm` instead of `plm`
pre <- lm(
  dflfp ~ flfplag + dflfplag +
    log_wdi_overvaluedlag + dlog_wdi_overvalued +
    log_gdpcaplag + dlog_gdpcap +
    sqlog_gdpcaplag + dsqlog_gdpcap +
    log_oilgascaplag + dlog_oilgascap +
    as_factor(irr_rate_regimelag) +
    wdi_fertilitylag + dwdi_fertility +
    vdem_elec_demlag + dvdem_elec_dem +
    as_factor(cowcode) - cowcode - year - country_space,
  data = df
)

# Take the data used in the estimation and put it in a data frame with tidy
# names
df2 <- pre$model %>% data.frame() %>% clean_names()

# Make the spatial weights where k = 10
W <- make_ntspmat(pre, ci = country_space, y = year, k = 10)

# Now, we can estimate the SAR model (may take a couple of minutes to run)
m_flfp_sar <- lagsarlm(
  dflfp ~ flfplag + dflfplag +
    log_wdi_overvaluedlag + dlog_wdi_overvalued +
    log_gdpcaplag + dlog_gdpcap +
    sqlog_gdpcaplag + dsqlog_gdpcap +
    log_oilgascaplag + dlog_oilgascap +
    as_factor_irr_rate_regimelag +
    wdi_fertilitylag + dwdi_fertility +
    vdem_elec_demlag + dvdem_elec_dem +
    as_factor(cowcode),
  listw = mat2listw(W[[2]], style = "W"),
  Durbin = FALSE,
  data = df2
)

# View the results
summary(m_flfp_sar)

# Table 2: FLFP ~ REER  + Country Effects    | REER @ 80% HDI  ----------

# What's the highest density 80% of the distribution?
# You may need to install the package `bayestestR` package for the code below

df_tscs %>% 
  filter(oecd == 0 & country != "Zimbabwe") |>
  summarize(hdi = bayestestR::hdi(log_wdi_overvalued, ci = .8))


m_flfp_strict <- plm(
  dflfp ~ flfplag + dflfplag +
    log_wdi_overvaluedlag + dlog_wdi_overvalued +
    log_gdpcaplag + dlog_gdpcap +
    sqlog_gdpcaplag + dsqlog_gdpcap +
    log_oilgascaplag + dlog_oilgascap - country,
  data = subset(
    df_tscs,
    oecd == 0 &
      country != "Zimbabwe" &
      log_wdi_overvalued >= -.48 &
      log_wdi_overvalued <= .64
  ),
  model = "within",
  effect = "individual",
  index = c("cowcode", "year")
)

# View the results
show_pcse(m_flfp_strict)

# Table 2: FLFP ~ REER  + Extended Controls  | REER @ 80% HDI  ----------

# First, re-create the relevant data
df <- df_tscs %>% 
  filter(
    oecd == 0 &
      country_space != "Zimbabwe (Rhodesia)" &
      country_space != "Tonga" &
      log_wdi_overvalued >= -.48 &
      log_wdi_overvalued <= .64
  ) 

# Estimate the model
pre <- lm(
  dflfp ~ flfplag + dflfplag +
    log_wdi_overvaluedlag + dlog_wdi_overvalued +
    log_gdpcaplag + dlog_gdpcap +
    sqlog_gdpcaplag + dsqlog_gdpcap +
    log_oilgascaplag + dlog_oilgascap +
    as_factor(irr_rate_regimelag) +
    wdi_fertilitylag + dwdi_fertility +
    vdem_elec_demlag + dvdem_elec_dem +
    as_factor(cowcode) -
    cowcode - country - country_space - year,
  data = df
)

# Create the usable data frame for the SLX model
df2 <- pre$model %>% data.frame() %>% clean_names()

# Make the spatial weights where k = 10
W <- make_ntspmat(pre, ci = country_space, y = year, k = 10)

# Estimate the SAR model (may take a couple of minutes)
m_flfp_strict_sar <- lagsarlm(
  dflfp ~ flfplag + dflfplag +
    log_wdi_overvaluedlag + dlog_wdi_overvalued +
    log_gdpcaplag + dlog_gdpcap +
    sqlog_gdpcaplag + dsqlog_gdpcap +
    log_oilgascaplag + dlog_oilgascap +
    as_factor_irr_rate_regimelag +
    wdi_fertilitylag + dwdi_fertility +
    vdem_elec_demlag + dvdem_elec_dem +
    as_factor(cowcode),
  listw = mat2listw(W[[2]], style = "W"),
  Durbin = FALSE,
  data = df2
)

# View the results
summary(m_flfp_strict_sar)


# Table 2: Ratio ~ REER + Country Effects ---------------------------------------------------

m_ratio_base <-
  plm(
    dflfp_mlfp ~ flfp_mlfplag + dflfp_mlfplag +
      log_wdi_overvaluedlag + dlog_wdi_overvalued +
      log_gdpcaplag + dlog_gdpcap + sqlog_gdpcaplag +
      dsqlog_gdpcap + log_oilgascaplag + dlog_oilgascap,
    data = subset(df_tscs, oecd == 0 & country != "Zimbabwe"),
    model = "within",
    effect = "individual",
    index = c("cowcode", "year")
  )

# View results with PCSEs
show_pcse(m_ratio_base)

# Table 2: Ratio ~ REER + Extended Controls ------------------------

# Create the usable dataset
df <- df_tscs %>%  
  filter(oecd == 0 &
           country_space != "Zimbabwe (Rhodesia)" &
           country_space != "Tonga")

pre <- lm(
  dflfp_mlfp ~ flfp_mlfplag + dflfp_mlfplag +
    log_wdi_overvaluedlag + dlog_wdi_overvalued +
    log_gdpcaplag + dlog_gdpcap +
    sqlog_gdpcaplag + dsqlog_gdpcap +
    log_oilgascaplag + dlog_oilgascap +
    as_factor(irr_rate_regimelag) +
    wdi_fertilitylag + dwdi_fertility +
    vdem_elec_demlag + dvdem_elec_dem +
    as_factor(cowcode) - cowcode - year - country_space,
  data = df
)

df2 <- pre$model %>% data.frame() %>% clean_names()

# Make the spatial weights where k = 10
W <- make_ntspmat(pre, ci = country_space, y = year, k = 10)

# Estimate the SAR model (may take a couple of minutes)
m_ratio_sar <- lagsarlm(
  dflfp_mlfp ~ flfp_mlfplag + dflfp_mlfplag +
    log_wdi_overvaluedlag + dlog_wdi_overvalued +
    log_gdpcaplag + dlog_gdpcap +
    sqlog_gdpcaplag + dsqlog_gdpcap +
    log_oilgascaplag + dlog_oilgascap +
    as_factor_irr_rate_regimelag +
    wdi_fertilitylag + dwdi_fertility +
    vdem_elec_demlag + dvdem_elec_dem +
    as_factor(cowcode),
  listw = mat2listw(W[[2]], style = "W"),
  Durbin = FALSE,
  data = df2
)

# View the results
summary(m_ratio_sar)


# Table 2: Ratio ~ REER + Country Effects   | REER @ 80% HDI  ----------

m_ratio_strict <- plm(
  dflfp_mlfp ~ flfp_mlfplag + dflfp_mlfplag +
    log_wdi_overvaluedlag + dlog_wdi_overvalued +
    log_gdpcaplag + dlog_gdpcap +
    sqlog_gdpcaplag + dsqlog_gdpcap +
    log_oilgascaplag + dlog_oilgascap,
  data = subset(
    df_tscs,
    oecd == 0 & country != "Zimbabwe" &
      log_wdi_overvalued >= -.48 &
      log_wdi_overvalued <= .64
  ),
  model = "within",
  effect = "individual",
  index = c("cowcode", "year")
)

# View the results
show_pcse(m_ratio_strict)


# Table 2: Ratio ~ REER + Extended Controls | REER @ 80% HDI  ----------

# Create the usable dataset
df <- df_tscs %>% 
  filter(
    oecd == 0 &
      country_space != "Zimbabwe (Rhodesia)" &
      country_space != "Tonga" &
      log_wdi_overvalued >= -.48 &
      log_wdi_overvalued <= .64
  )

# Estimate the preliminary model
pre <- lm(
  dflfp_mlfp ~ flfp_mlfplag + dflfp_mlfplag +  
    log_wdi_overvaluedlag + dlog_wdi_overvalued +
    log_gdpcaplag + dlog_gdpcap +
    sqlog_gdpcaplag + dsqlog_gdpcap +
    log_oilgascaplag + dlog_oilgascap +
    as_factor(irr_rate_regimelag) + 
    wdi_fertilitylag + dwdi_fertility +
    vdem_elec_demlag + dvdem_elec_dem +
    as_factor(cowcode) - cowcode - year - country_space,
  data = df
)

# Make the data frame suitable for the SLX model
df2 <- pre$model %>% data.frame() %>% clean_names()

# Make the spatial weights where k = 10
W <- make_ntspmat(pre, ci = country_space, y = year, k = 10)

# Estimate the SAR model (may take a couple of minutes)
m_ratio_strict_sar <- lagsarlm(
  dflfp_mlfp ~ 0 + flfp_mlfplag + dflfp_mlfplag + 
    log_wdi_overvaluedlag + dlog_wdi_overvalued +
    log_gdpcaplag + dlog_gdpcap +
    sqlog_gdpcaplag + dsqlog_gdpcap +
    log_oilgascaplag + dlog_oilgascap +
    as_factor_irr_rate_regimelag + 
    wdi_fertilitylag + dwdi_fertility +
    vdem_elec_demlag + dvdem_elec_dem +
    as_factor(cowcode),
  listw = mat2listw(W[[2]], style = "W"),
  Durbin = FALSE,
  data = df2
)

# View the results
summary(m_ratio_strict_sar)

# Table 2: Make the table ----------------------------------------------------

# Sanity test. Do my custom extract functions extract the right standard errors
# and p-values
coeftest(m_flfp_base, vcov. = vcovBK(m_flfp_base, type = "HC1", cluster = "time"))
show_pcse(m_flfp_base)
data.frame("pcse" = extract_pcse(m_flfp_base), 
           "pcse_sig" = extract_pcse_sig(m_flfp_base))

screenreg(
  list(
    m_flfp_base,
    m_flfp_sar,
    m_flfp_strict,
    m_flfp_strict_sar,
    m_ratio_base,
    m_ratio_sar,
    m_ratio_strict,
    m_ratio_strict_sar
  ),
  override.se = list(
    extract_pcse(m_flfp_base),
    m_flfp_sar,
    extract_pcse(m_flfp_strict),
    m_flfp_strict_sar,
    extract_pcse(m_ratio_base),
    m_ratio_sar,
    extract_pcse(m_ratio_strict),
    m_ratio_strict_sar
  ),
  override.pvalues = list(
    extract_pcse_sig(m_flfp_base),
    m_flfp_sar,
    extract_pcse_sig(m_flfp_strict),
    m_flfp_strict_sar,
    extract_pcse_sig(m_ratio_base),
    m_ratio_sar,
    extract_pcse_sig(m_ratio_strict),
    m_ratio_strict_sar
  ),
  # Add significance stars
  stars = c(0.01, 0.05, 0.1),
  # Add a header 
  custom.header = list(
    "$\\Delta$FLFP" = 1:4,
    "$\\Delta$Labor Force Ratio" = 5:8
  ),
  # Variable names
  custom.coef.map = list_var_names_texreg,
  # Goodness-of-fit statistics
  include.rsquared = FALSE,
  include.adj = TRUE,
  include.adjrs = FALSE,
  include.aic = FALSE,
  include.loglik = FALSE,
  include.lr = FALSE,
  # Custom notes
  custom.note = "%stars,
  All models include country-specific fixed effects",
)



# LaTeX Code
# texreg(
#   list(
#     m_flfp_base,
#     m_flfp_sar,
#     m_flfp_strict,
#     m_flfp_strict_sar,
#     m_ratio_base,
#     m_ratio_sar,
#     m_ratio_strict,
#     m_ratio_strict_sar
#   ),
#   override.se = list(
#     extract_pcse(m_flfp_base),
#     m_flfp_sar,
#     extract_pcse(m_flfp_strict),
#     m_flfp_strict_sar,
#     extract_pcse(m_ratio_base),
#     m_ratio_sar,
#     extract_pcse(m_ratio_strict),
#     m_ratio_strict_sar
#   ),
#   override.pvalues = list(
#     extract_pcse_sig(m_flfp_base),
#     m_flfp_sar,
#     extract_pcse_sig(m_flfp_strict),
#     m_flfp_strict_sar,
#     extract_pcse_sig(m_ratio_base),
#     m_ratio_sar,
#     extract_pcse_sig(m_ratio_strict),
#     m_ratio_strict_sar
#   ),
#   # Add significance stars
#   stars = c(0.01, 0.05, 0.1),
#   # Add a header 
#   custom.header = list(
#     "$\\Delta$FLFP" = 1:4,
#     "$\\Delta$Labor Force Ratio" = 5:8
#   ),
#   # Variable names
#   custom.coef.map = list_var_names_texreg,
#   # Goodness-of-fit statistics
#   include.rsquared = FALSE,
#   include.adj = TRUE,
#   include.adjrs = FALSE,
#   include.aic = FALSE,
#   include.loglik = FALSE,
#   include.lr = FALSE,
#   # Custom notes
#   custom.note = "%stars,
#   All models include country-specific fixed effects",
#   # LaTeX formatting options
#   booktabs = TRUE,
#   fontsize = "scriptsize",
#   siunitx = TRUE,
#   threeparttable = FALSE,
#   sideways = TRUE,
#   float.pos = "tp",
#   label = "tab:results"
# )

#################################################################
##  Table 3: Female Employment ~ Industry | country == Mexico  ##
#################################################################


# Table 3: Make data: ILO | Mexico -----------------------------

# First, make a dataset of ILO variables just for the case of Mexico. 
# Also, create logged and lagged values
df_ilo_ind_mex <- df_ilo_ind %>% 
  filter(country == "Mexico") %>% 
  arrange(industry, year) %>% 
  mutate(log_employment = log(employment)) %>% 
  group_by(industry) %>% 
  mutate(log_employmentlag = dplyr::lag(log_employment, 1)) %>% 
  ungroup()
# Table 3: Female Employment ~ ILO Industry Fixed Effects | country == "Mexico" ---------

# Make a panel object for use with `feols`
pdat <- panel(df_ilo_ind_mex, ~ industry + year)

# Estimate the model. Note: Manufacturing is the excluded category
m_emp_ind_mex <- feols(log_employment ~ industry, data = pdat)

# Footnote 10 refers to various robustness checks. Estimate them here.
m_emp_ind_mex_footnote10_v1 <-
  feols(log_employment ~ industry + l(log_employment, 1:2), data = pdat)

m_emp_ind_mex_footnote10_v2 <-
  feols(log_employment ~ industry + year + l(log_employment, 1:2), data = pdat)

m_emp_ind_mex_footnote10_v3 <-
  feols(log_employment ~ industry + factor(year) + l(log_employment, 1:2),
        data = pdat)

# View the results with clustered standard errors
summary(m_emp_ind_mex, cluster = c("industry", "year"))
summary(m_emp_ind_mex_footnote10_v1, cluster = c("industry", "year"))
summary(m_emp_ind_mex_footnote10_v1, cluster = c("industry", "year"))
summary(m_emp_ind_mex_footnote10_v1, cluster = c("industry", "year"))

# Table 3: Make the table  -------------------------------------------------

# Create a list of variable names for the texreg table
list_ilo_var_names_texreg <- list(
  'industryAccommodation and food service'         = 'Accommodation and food service',
  'industryAdministrative and support services'    = 'Administrative and support services',
  'industryAgriculture; forestry and fishing'      = 'Agriculture; forestry and fishing',
  'industryArts, entertainment and recreation'     = 'Arts, entertainment and recreation', 
  'industryConstruction'                           = 'Construction',  
  'industryEducation'                              = 'Education',  
  'industryElectricity; gas, and steam'            = 'Electricity; gas, and steam',  
  'industryExtraterritorial organizations'         = 'Extraterritorial organizations',
  'industryFinance and insurance'                  = 'Finance and insurance',  
  'industryHealth and social work'                 = 'Health and social work',  
  'industryActivities of household as employers'   = 'Household as employers',
  'industryInformation and communication'          = 'Information and communication',  
  'industryMining and quarrying'                   = 'Mining and quarrying',
  'industryNot elsewhere classified'               = 'Not elsewhere classified',
  'industryOther services'                         = 'Other services',
  'industryPublic administration'                  = 'Public administration', 
  'industryReal estate activities'                 = 'Real estate activities',
  'industryScientific and technical activities'    = 'Scientific and technical activities',
  'industryTransportation and storage'             = 'Transportation and storage',
  'industryWater supply; sewage, waste management' = 'Water supply; sewage, waste management',
  "industryWholesale and retail trade"             = 'Wholesale and retail trade',
  'year'                                           = 'Time Trend', 
  'log_employmentlag'                              = 'Female Employment (Logged)$_{(t-1)}$',
  '(Intercept)'                                    = 'Constant'
)

# make the table
screenreg(
  m_emp_ind_mex,
  # Get the standard errors and p-values right
  override.se = summary(m_emp_ind_mex, cluster = c("industry", "year"))$coeftable[, 2],
  override.pvalues = summary(m_emp_ind_mex, cluster = c("industry", "year"))$coeftable[, 4],
  single.row = TRUE,
  stars = c(0.01, 0.05, 0.1),
  caption = "",
  # Coefficient labels
  custom.coef.map = list_ilo_var_names_texreg,
  # Goodness-of-fit statistics
  include.rsquared = FALSE,
  include.adjrs = FALSE,
  include.loglik = FALSE,
  # Custom notes
  custom.note = "%stars,
  Standard errors clustered by industry and year in parentheses",
)

#################################################################
##               Table 4: FLFP Models for Mexico               ##
#################################################################

# Table 4: Make data: FLFP ~ REER | Mexico ----------------------------------------

# Create a database of that just the case of Mexico since 1990.
df_mex_case <- df_tscs %>% filter(country == "Mexico" & year >= 1990)

# Table 4: FLFP ~ REER | Mexico ---------------------------------------------

# Estimate the ECM model via `dynardl`
m_mex_flfp <- dynardl(
  flfp ~ log_wdi_overvalued + log_gdpcap,
  data = df_mex_case,
  lags = list(
    "flfp" = 1,
    "log_wdi_overvalued" = 1,
    "log_gdpcap" = 1
  ),
  diffs = c("log_wdi_overvalued",
            "log_gdpcap"),
  lagdiffs = list("flfp" = 1:2),
  trend = FALSE,
  ec = TRUE,
  simulate = FALSE, # No need to simulate here
  fullsims = FALSE,
  shockvar = "log_wdi_overvalued",
  shockval = -(sd(df_mex_case$log_wdi_overvalued)),
  time = 2,
  range = 10
)

# View the results
summary(m_mex_flfp)

# Table 4: Labor Force Ratio ~ REER | Mexico ---------------------------------------------

# Estimate the ECM model via `dynardl`
m_mex_flfp_mlfp <- dynardl(
  flfp_mlfp ~ log_wdi_overvalued + log_gdpcap,
  data = df_mex_case,
  lags = list(
    "flfp_mlfp" = 1,
    "log_wdi_overvalued" = 1,
    "log_gdpcap" = 1
  ),
  diffs = c("log_wdi_overvalued",
            "log_gdpcap"),
  lagdiffs = list("flfp_mlfp" = 1:2),
  trend = FALSE,
  ec = TRUE,
  simulate = FALSE,
  fullsims = FALSE,
  shockvar = "log_wdi_overvalued",
  shockval = -(sd(df_mex_case$log_wdi_overvalued)),
  time = 2,
  range = 10
)

# View the results
summary(m_mex_flfp_mlfp)

# Table 4: Make the table ----------------------------------------------------------

# Make variable names
list_var_names_mex_texreg <- list(
  'l.1.flfp'                =  'LDV$_{(t-1)}$',
  'ld.1.flfp'               =  '$\\Delta$LDV$_{(t-1)}$',
  'ld.2.flfp'               =  '$\\Delta$LDV$_{(t-2)}$',
  'l.1.flfp_mlfp'           =  'LDV$_{(t-1)}$',
  'ld.1.flfp_mlfp'          =  '$\\Delta$LDV$_{(t-1)}$',
  'ld.2.flfp_mlfp'          =  '$\\Delta$LDV$_{(t-2)}$',
  'l.1.log_wdi_overvalued'  =  'logOvervalued$_{(t-1)}$',
  'd.1.log_wdi_overvalued'  =  '$\\Delta$logOvervalued$_t$', 
  'l.1.log_gdpcap'          =  'GDP$_{(t-1)}$',
  'd.1.log_gdpcap'          =  '$\\Delta$GDP$_t$',
  '(Intercept)'             = 'Constant'
)

# make the table
screenreg(
  list(m_mex_flfp,
       m_mex_flfp_mlfp),
  stars = c(0.01, 0.05, 0.1),
  # Coefficient labels
  custom.coef.map = list_var_names_mex_texreg,
  # Custom notes
  custom.note = "%stars",
)


#################################################################
##                 Table 5: Political Interest                 ##
#################################################################


# Table 5: Political Interest ---------------------------------------------

source("04_political_interest_results.R")

# View the table
tab_interest

#################################################################
##                           Figures                           ##
#################################################################

## Figure 1: Draw in Latex
# Figure 2: Overvalued Histogram --------------------------------------------

fig_hist_overvalued <-
  ggplot(data = df_tscs %>% filter(oecd == 0 &
                                     country != "Zimbabwe"),
         aes(x = log_wdi_overvalued)) +
  geom_histogram(bins = 100,
                 fill = "grey70",
                 color = "grey70") +
  labs(x = "logOvervalued", y = "Count")

fig_hist_overvalued


# Figure 3: Marginal Effect: logOvervalued ~ Interest * Regime Type -------

# Note, results require the `lm_interest` object created in the 

# Make the plot
fig_me_interest_dem <-
  plot_slopes(
    lm_interest,
    variables = "some_interest_share",
    condition  = "vdem_elec_dem",
    vcov = "HC1",
    rug = TRUE
  ) + 
  geom_hline(yintercept = 0, linetype = "dashed") + 
  labs(x = "Regime Type", 
       y = "Marginal EFfect \nof Political Interest")

fig_me_interest_dem


#################################################################
##                    Supplemental Appendix                    ##
#################################################################



# List of Countries -------------------------------------------------------

list_countries <- df_sum_stats %>%
  select(country) %>%
  group_by(country) %>% 
  slice(1) %>% 
  arrange(country) %>% 
  ungroup()

# View the country list 
data.frame(list_countries)

# LaTeX code was generated by hand

# SI: GMM Results  --------------------------------------------------------

cat(paste(
  "I conduct these tests in Stata using the following script:
  'scripts/03_stata_gmm_results.do' 
  Use that script to replicate the results."))



# Session Information -----------------------------------------------------

sessionInfo()
